Increased moist heat stress risk across China under warming climate

Heatwaves have afflicted human health, ecosystem, and socioeconomy and are expected to intensify under warming climate. However, few efforts have been directed to moist heat stress (MHS) considering relative humidity and wind speed, and moist heat stress risk (MHSR) considering exposure and vulnerability. Here we showed MHS and MHSR variations across China during 1998–2100 using China Meteorological Administration Land Data Assimilation System datasets, the 6th Coupled Model Intercomparison Project (CMIP6) merged datasets, Gross Domestic Product, population and leaf area index. We detected increased MHS across China under different Shared Socioeconomic Pathways (SSPs). Specifically, the historical MHS occurred mostly during mid-July to mid-August. We found increasing trends of 0.08%/year, 0.249%/year, and 0.669%/year in the MHS-affected areas under SSP126, SSP245, and SSP585, respectively. Furthermore, we observed the highest increasing rate of MHSR in Northwest and Southwest China, while the MHSR across Northeast and North China under SSP126 shifted from increasing to decreasing trends. Noteworthy is that the increasing trend of MHSR under SSP585 is 1.5–2.6 times larger than that under SSP245, especially in North and South China. This study highlights spatiotemporal evolutions of MHS and MHSR and mitigation to moisture heat stress in a warming climate.

www.nature.com/scientificreports/ of MHSR across China can greatly help to enhance mitigation to MHS at regional scales. This point justify the significance of this current study. The objectives of this study therefore are to: (1) realize the data fusion by fully considering the systematic error of the CMIP6 model to improve its modelling accuracy; (2) delineate spatiotemporal evolutions of MHS considering the influences of RH and WS over China; and (3) quantify the spatiotemporal evolutions of regional MHSR across China considering population, GDP and leaf area index (LAI). The study seeks to foster improved understanding of MHS evolution and enhancement of mitigation of MHSR in a warming climate.

Results and discussion
Evaluation of CLDAS and CMIP6 datasets. To ensure calculation accuracy of humidity (HI) and Steadman Index (SI), we evaluated the quality of the CLDAS data using the observed data of CMA (The distribution of observation stations can be seen in Fig. S1) and compared the CLDAS with the Global Land Data Assimilation System dataset (GLDAS) 28 and the fifth-generation ECMWF reanalysis datasets (ERA5) 29 The bias, RMSE, and Corr. of ERA5, GLDAS, and CLDAS with in-situ observed daily maximum temperature are shown in Fig. 1. In general, the bias of CLDAS daily maximum temperature was smaller than that of ERA5 and GLDAS daily maximum temperature. The RMSE of CLDAS was 0.9 °C, which was distinctly smaller than ERA5 (~ 2.5 °C) and GLDAS (~ 2.9 °C), and the Corr. between CLDAS temperature data and the in-situ observed temperature was about 0.98, being significantly larger than ERA5 (~ 0.9) and GLDAS (~ 0.86). The evaluation results for the specific humidity and wind speed can be found in Fig. S5 and Fig. S6 of the Supplementary Information. Therefore, the CLDAS data performed better than ERA5 and GLDAS data in describing changes in the summer daily maximum temperature and the specific humidity across China. Yin et al. used data such as air temperature fused with station observations to analyze apparent temperature in the one belt and one road region, whose the RMSE of air temperature is between 1 and 2 °C 16 , which indirectly verified the accuracy of our data source.
To ensure the accuracy of the CMIP6 data after the spatial downscaling and data fusion. Here we calculated the bias between CLDAS and the original FGOALS and CanESM outputs, the SSD-based processed and the   (Fig. 2). We found a large positive bias (~ 2-6 °C) between the original FGOALS data and the CLDAS in Eastern China, and a large negative bias (− 6 to − 2 °C) in West China (Fig. 2a). We aslo detected a negative bias (− 3 to − 1 °C) between CanESM outputs and CLDAS in South China (Fig. 2b), a positive bias (2-3 °C) in most parts of North and Northeast China, a larger positive bias in the middle of the Tibetan Plateau (4 to 6 °C), and a large negative bias in the western part of the Tibetan Plateau (− 6 to − 2 °C). The bias of the FGOALS and CanESM was reduced using the SSD method, and the bias was − 1 to 1 °C (Fig. 2c,d). The TC-based on fusion determines the weight of the models based on the covariance of each model. Here we also included advantages of FGOALS and CanESM data after SSD-based data processing. The bias was small in some regions such as North China and Northeast China. The bias of relative humidity can be found in Fig. S7 of the Supplementary Information, and the effect was the same as for air temperature, which was better than the original data by SSD-based and TC-based data processing.    The spatial evolution of the SI MHS intensity considering wind speed indicates that the changing magnitude of SI MHS intensity was lower than HI MHS intensity (Fig. S11). Under SSP126 and SSP245 scenarios, the SI heatwave intensity in the Tibetan Plateau, North China, and Northeast China evidently increased and then followed by moderate changes. The SI heatwave intensity in China would generally increase under the SSP585 scenario and the largest increasing magnitude of SI heatwave intensity was found in the Tibetan Plateau, North China, and Northeast China. Little changes in SI heatwave intensity were observed in South China.
These results indicated increased frequency ( Fig. S9 and S10), intensity (   We integrated the exposure of MHS, population, GDP and LAI under different scenarios to represent the spatiotemporal evolutions of MHSR. Yin et al. analyzed the spatiotemporal distribution and risk of heat waves by apparent temperature, population and NDVI data, and the results showed that the eastern part of China deepens economic development while taking heat wave risk into account 31 . Here we better reflected the spatial variations by CEEMDAN method under different scenarios. While GDP is also an important factor in heat stress risk for low-income people. There are relatively few countermeasures for heat stress, for example, expensive air conditioners are a burden for them 32 , and the level of economic development in North China is relatively low compared to Southeast China. Therefore it can be seen that the spatial extent of MHSR in North China is wider compared to Southeast China. To further analyze the temporal trends of MHSR, we analyzed the temporal changes of MHSR in different regions in China using the MK test. For Northeast China (Fig. 6a), MHSR first increased and then showed a decreasing trend under SSP126, with an overall increase of 0.00049/year, while MHSR showed a continuously increasing trend under SSP245 and SSP585, with 0.00108/year and 0.00308/year, respectively. For North China www.nature.com/scientificreports/ (Fig. 6b), the MHSR slowly increased under SSP126, with an overall increase of 0.00049/year, and the MHSR also showed a continuous increase under SSP245 and SSP585, with 0.00108/year and 0.00308/year, respectively. For South China (Fig. 6c), the MHSR first increased and then stabilized under SSP126, with an overall increase of 0.00037/year, while the MHSR showed a continuous increase under SSP245 and SSP585, with 0.00233/year and 0.00385/year, respectively. For Northwest China (Fig. 6d), the MHSR increased and then decreased under SSP126, with an overall increase of 0.00145/year, while the MHSR showed a continuous increase under SSP245 and SSP585, with 0.00308/year and 0.00462/year, respectively. For Southwest China (Fig. 6e), MHSR increased and then decreased under SSP126, with an overall increase of 0.001408/year, while MHSR showed a continuous increasing trend under SSP245 and SSP585, with 0.00282/year and 0.00423/year, respectively. For the Tibetan Plateau (Fig. 6f)

Discussion
In the context of global warming, the burning of fossil fuels, deforestation, rapid urbanization, growing population brought about by economic development and the anthropogenic heat generated by the overuse of energy, people will face more and more MHS in the future, especially for China, a developing country with a huge population 10,11 . The frequency of moist heat stress events, increased population exposure, and increasing human health problems caused by MHS have affected economic development to some extent, especially for urban residents with poor ventilation, buildings and impervious surfaces that absorb and retain heat, as well as increased high humidity under anthropogenic heat release 16,24,25 .
Nowadays, more and more countries are taking corresponding measures to reduce the MHSR, such as urban greening, irrigation, central air conditioning and the construction of corresponding heat avoidance facilities 7,21,27 . Urban greening increases the leaf area index, absorbing more solar radiation and reducing the temperature increase due to absorption of solar radiation on the ground. While the transpiration of vegetation can effectively improve the heat-sensing and latent heat changes in the subsurface, regulating the relative humidity and reducing heat stress to a certain extent. Development of socioeconomy can allow air conditioning and which can greatly alleviate indoor heat stress-induced healthy issues. While, energy consumption also generates myriad anthropogenic heat 22,23,27 . Here we found changing pattern of MHSR as "increases-decreases" under SSP126 scenarios in most regions, while for SSP245 and SSP585 scenarios, the increasing trend of MHSR under SSP585 scenarios is 1.5-2.6 times larger than that of SSP245 scenarios. So the sustainable development path can help to reduce the MHSR such as the carbon neutrality and carbon peaking policies currently being implemented in China.
In terms of data source uncertainty, the CLDAS air temperature, relative humidity and wind speed fused with more than 50,000 in situ observations from the China Meteorological Administration (CMA) have a higher accuracy than ERA5 and GLDAS. However, sparse distribution of CMA station observations in western China when compared to the eastern China results in relatively lower quality of CLDAS data in the western than in the eastern China 33 . Besides, we spatially downscaled and fused the CMIP6 data, uncertainty is still unavoidable in the CMIP6 data under different climate models due to different model types, parameterization schemes, and spatial resolutions 14,31 . What's more, the GLASS leaf area index data we used from 1998-2018 are mainly based on the results of satellite inversion, which have relatively higher accuracy and are better than MODIS and AVHRR. Certain errors can also be expected for areas with remarkable heterogeneous vegetation cover 34,35 . The CMIP6 data during 2020-2100 used for the LAI data under different scenarios are mainly generated by climate models, which have higher quality than the CMIP5 model, but there is an overestimation of LAI in non-forested regions 36 .

Conclusions
In this study, we depicted MHS and MHSR changes in both space and time across China during a period of 1998-2100 based on CLDAS datasets, CMIP6 datasets, GDP, population and LAI datasets under SSP126, SSP245, and SSP585 scenarios. We have obtained the following important and interesting conclusions: Considering the coarse spatial resolution and systematic errors of the models in CMIP6, we adopt the SD and TC methods to realize the spatial downscaling of the two models and the calculation of systematic errors for each grid point of the models, respectively. Then we calculate their weights according to the model errors and realize the fusion of CMIP6 model data, so as to reduce the uncertainty brought by the data source and ensure the accuracy and quality of the data source. It can be found that the fusion results are significantly better than the original data, especially the bias was obviously reduced.
The MHS in China showed an increasing temporal trend in a warmer climate. Firstly, the MHS-affected area in China has gradually increased, with an average annual increase of 0.399% and 0.479% under HI and SI from 1998 to 2020. The MHS-affected area under SSP126 first increased and then stabilized from 2021 to 2100. The MHS under SSP245 and SSP585 were gradually increasing, and the MHS-affected area under SSP585 was higher than that under SSP126 and SSP245. Besides, the frequency of historical MHS mainly occurred from mid-July to mid-August. The MHS under SSP126 also mainly occurred from early July to early August, while the MHS under SSP245 and SSP585 would continue from early July to the end of August. Meanwhile, the frequency under SSP585 was higher than that under SSP126 and SSP245.
The MHSR in China showed an increasing temporal trend in a warmer climate except for Northeast China and North China under SSP126. Firstly, it can be seen from the spatial evolution that the extent and magnitude All in all, we analyzed the intensity, affected area, duration, and frequency of MHS in China, and explored the spatial and temporal evolution of MHSR in different regions of China considering the hazards, exposures, and vulnerability by combining population, GDP, and LAI data. These results can provide a basis for the management and response to moist heat stress risk in different regions.

Method
The spatial downscaling and data fusion of the CMIP6 datasets. To ensure the accuracy of the CMIP6 data used, we used the spatial disaggregation (SD) statistical downscaling method for spatial downscaling of the CMIP6 data 37 . Considering that the bias in monthly average temperature by the traditional SD method would smooth information, especially the temperature, we improved a sliding SD method (SSD) at a three-day scale (Fig.S2). We also fully considered and calculated the error of the model using the Triple Collocation method (TC) 38 . Then we calculated the weights of model according to the error of model and obtained the fusion of FGOALS and CanESM models.
The variance and covariance formulas for CLDAS, FGOALS and CanESM models are shown as: The systematic error of the CLDAS, FGOALS and CanESM models data sets are as follows: The weights of each model were obtained according to the error of each model, and the weights were calculated as: Multi-model fusion was performed according to the weights of different models of CMIP6, and the CMIP6 results after multi-model fusion can be obtained: The spatial downscaling and model-fusion procedure was shown in Fig. S2 and the sample figure was in Fig. S3. The calculation of MHSR. According to the risk triangle evaluation theory, the calculation of MHSR mainly takes into account hazard, exposure and vulnerability 21 . Hazard mainly considers the frequency, duration and intensity of MHS, exposure mainly considers the impact of MHS on population and GDP, and vulnerability mainly considers the mitigation effect of higher GDP and vegetation cover on MHS. All data were normalized before the MHSR was calculated, and the MHSR equation is as follows. where T denotes the air temperature, RH denotes the relative humidity, P is vapor pressure, and V denotes the wind speed.
Here we used relative and absolute temperature thresholds to identify MHS and 90% percentile was taken as the relative temperature threshold for MHS 13,41 . The MHS of 3 °C higher than the long-term mean MHS is taken as the absolute temperature threshold. Given that the relative MHS threshold is lower than the absolute MHS threshold, the absolute MHS threshold can be accepted to define heatwaves 16 .The spatial distribution of the thresholds of the MHS can be found in Fig. S4 in Supplementary Information. Spatial-temporally varying trend. The Mann-Kendall (MK) test was used to detect the trends in the meteorological serie 42 . The separate spatiotemporally-varying trend of the heatwave from 1998 to 2100 was diagnosed using the Complete Ensemble Empirical Mode Decomposition with Adaptive Noise Analysis method (CEEMDAN) 43 . We diagnosed the value increment of the CEEDMAN spatial trend at a given time from 1998 to 2100, and the reference time was 1998 year, which represented accumulated warming from 1998: H(t) is the original signal, IMF k is the Intrinsic Mode Function, and R(t) is the extracted trend.

Data availability
The China Meteorological Administration Land Data Assimilation System (CLDAS) air temperature, specific humidity, and wind data with daily time scale and 0.25° spatial resolution from 1998 to 2020 are freely available at http:// data. cma. cn 44 . The Coupled Model Intercomparison Project datasets (CMIP6) air temperature, specific humidity, wind data and LAI data under SSP126, SSP245, and SSP585 scenarios from 2021 to 2100 are freely available at https:// esgf-node. llnl. gov/ search/ cmip6/ 45 . The population and GDP from 2000 to 2019 are freely available at https:// www. resdc. cn. The population and GDP data from 2020 to 2100 under SSP126, SSP245, and SSP585 scenarios are freely available at https:// geogr aphy. nuist. edu. cn 46 LAI datasets used the Global Land Surface Satellite (GLASS) LAI product from 1998 to 2018 are freely available at http:// glass-produ ct. bnu. edu. cn/ intro ducti on/ LAI. html 47 . The datasets used and analysed during the current study available from the corresponding author on reasonable request.